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Abstract 

One crucial feature of zygotic linkage disequilibrium (LD) analysis is its direct use of diploid genotyping data, irrespective of 
the type of mating system. Previous theories from an evolutionary perspective mainly focus on gametic LD, but the 
equivalent development for zygotic LD is not available. Here I study the evolution of zygotic LD and the covariances 
between gametic and zygotic LDs or between distinct zygotic LDs in a finite local population under constant immigration 
from a continent population. I derive the analytical theory under genetic hitchhiking effects or in a neutral process. Results 
indicate that zygotic LDs (diploid level) are more informative than gametic LD (haploid level) in indicating the effects of 
different evolutionary forces. Zygotic LDs may be greater than or comparable to gametic LD under the epistatic selection 
process, but smaller than gametic LD under the non epistatic selection process. The covariances between gametic and 
zygotic LDs are strongly affected by the mating system, linkage distance, and genetic drift effects, but weakly affected by 
seed and pollen flow and natural selection. The covariances between different zygotic LDs are generally robust to the 
effects of gene flow, selection, and linkage distance, but sensitive to the effects of genetic drift and mating system. 
Consistent patterns exist for the covariances between the zygotic LDs for the two-locus genotypes with one common 
genotype at one locus or without any common genotype at each locus. The results highlight that zygotic LDs can be 
applied to detecting natural population history. 
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Introduction 

Zygotic linkage disequilibrium (LD) refers to the difference 
between the joint genotypic frequency at two loci and the product 
of genotypic frequencies at each locus [1], [2], [3], [4], [5], [6]. 
The concept itself is a purely statistical term, and can also be 
viewed as the covariance of genotypic frequencies, analogous to 
the covariance of allelic frequencies for the concept of gametic LD 
[7] , [8] , [9] . Its biological significance can be viewed when used for 
detecting the effects of evolutionary forces by comparing its 
empirical distribution with the predicted distribution once an 
evolutionary model is specified [9], [10]. The commonality 
between gametic and zygotic LDs lies in their utility for measuring 
non-random associations between loci. The crucial difference is 
that zygotic LD analysis does not require a random mating 
assumption since it directly uses diploid genotyping data. 
However, gametic LD calculation inferred from diploid genotypes 
needs this assumption since haplotypes must be priorly known. 
Such a difference is significant because the potential false-positive 
errors could be substantial in inferring haplotypes/linkage phases 
using the diploid genotyping data sampled from a natural 
population of a mixed mating system. 

Previously relevant theories emphasize the joint frequency of 
double heterozygotes or double homozygotes in a neutral process, 
or the joint descent measures for a population with a mixed mating 
system [3], [4]. [1 1], [12], [13], [14], [15]. Zygotic LD is implicitly 
indicated from interpreting character associations in a partial 
inbreeding system [16], or from explaining an excess of the 



equilibrium genotypic frequencies at two independent loci in a 
mixed mating system [17], or from defining the covariance of 
heterozygosities [18] or the covariance of descent identities [19]. 
More recent studies concentrate on the statistical issues, including 
the procedure for testing zygotic LDs [5], [6], [20], [21], [22] and 
the potential application of zygotic LD to mapping quantitative 
trait loci (QTL) [23]. Unlike gametic LD that has received 
considerably theoretical studies from the evolutionary perspective 
[9], [10], [24], an equivalent theory for zygotic LD has not been 
fully developed. Although the evolutionary forces acting on the 
gametic LD may, in principle, also affect the zygotic LD, these 
effects and the resultant patterns have not been explicitly studied. 
This void motivates me to study how zygotic and other high-order 
LDs evolve under the effects of different evolutionary forces. 

In flowering plants, three distinct processes in a life cycle are 
involved in changing zygotic LD and its relationship with gametic 
LD in a local population. One process is the asymmetric 
immigration through haploid pollen flow and diploid seed flow 
from a source population. Pollen flow directly generates gametic 
LD, but indirectly affects zygotic LD since each pollen grain only 
carries one gamete in fusion with ovules in the recipient 
population. Seed flow can generate both zygotic and gametic 
LDs since each seed carries two gametes into the recipient 
population simultaneously. 

The second process influencing zygotic LD in plants is the 
mating system [25]. Selfing facilitates both gametic and zygotic 
LDs, even for the loci with a free recombination [17]; while 
random mating erodes both gametic and zygotic LDs. This effect 
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can be unequal on zygotic and gametic LDs since zygotic LD 
might be more sensitive than gametic LD to the mating system. 

The third process influencing zygotic LD in plants is selection in 
either the gametophyte or the sporophyte stage, or in both stages. 
Selection against heterozygote or epistatic selection at the 
sporophyte stage can directly change zygotic LD, but indirecdy 
changes gametic LD [26], such as in natural hybrid zones [27], 
[28], [29], [30]. This is analogous to the conventional artificial 
selection that directly exerts effects on zygotic LD but indirecdy on 
gametic LD in plant and animal breeding programs. Selection in 
the gametophyte stage direcdy changes gametic LD, but indirecdy 
changes zygotic LD owing to the connection between the 
gametophyte and sporophyte stages in one life cycle. The natural 
overloading of pollen on the stigma of a flower implies pollen 
competition and the occurrence of natural selection [31], [32], 
[33]. Also, an excess of ovule abortion in many single embryo or 
polyembryony plants implies the occurrence of natural selection in 
ovules [34], [35]. Some genes can express at both the gametophyte 
and sporophyte stages [36], [37], [38], [39], [40], and might 
experience different extents of selection pressure. Selection can 
change both gametic and zygotic LDs among these genes. Thus, it 
is of both theoretical and practical significance to study how the 
above three distinct processes evolve zygotic LD. 

Here I examine how different driving forces (mating system, 
genetic drift, migration, and natural selection) affect zygotic LD 
from the evolutionary perspective, complementary to the existing 
statistical issues. An island-continent model is considered, with an 
emphasis on the evolution of zygotic LD in the finite island 
population. I begin by presenting an exact model and use Monte 
Carlo (MC) simulations to evaluate the evolution of zygotic LD 
under different evolutionary forces in the population with a mixed 
mating system. I then derive analytical theories in two specific 
cases (genetic hitchhiking effects and a neutral process) under 
random mating, and validate the theories through MC simula- 
tions. Through the analytical and simulation results, I explore the 
evolutionary properties of zygotic LDs and discuss their potential 
utility. 

Results 

Exact Model 

Consider an island population, with constant immigration rates 
of pollen imp ) and seeds (mj) per generation from a continent 
population. The continent population is assumed to be sufficiendy 
large in size and stable in genetic composition. Migration from the 
island to continent population is neglected, and mutation effects 
are excluded. At the gametophyte stage, pollen and ovules are 
subject to natural selection before they are combined to produce 
seeds. The plant life cycle follows a sequence of events: pollen and 
ovules generation, pollen flow, selection at the gametophyte stage, 
mixed mating with a selfing rate a (0<a< 1), seed flow, selection 
at the sporophyte stage, and genetic drift. Selection strength may 
be either strong or weak; or epistatic selection is allowed in either 
the gametophyte or the sporophyte stage, or in both. The same 
mating system is assumed in the island and continent populations. 

Consider two diallelic nuclear loci, with alleles ^4iand A2 at 
locus A, alleles B\ and B2 at locus B, and a recombination rate r 
between the two loci. These alleles may refer to single nucleotide 
polymorphic (SNP) markers since tri- or tetra-allelic SNP sites are 
infrequent in natural populations. In the island population, let 
PA,B t {i, k= 1, 2) and D ^#be the frequency of gamete AjB^ and the 
gametic LD in current adults, respectively; and PAiB k C3n be 
expressed as PAiB k =PAiPB k +(—^)' +k Dab- For the random 
mating part, the frequency of gamete AjBk in pollen or ovules 



(the next generation) produced by the current adults can 
be expressed as p* AiBi =p A .p Bl: +(- l)' +k (l -r)D AB [41]. 
Similarly, let PAiAjB k B,> Pa,Aj, and PB k B, be the frequencies 
of genotypes AjAjB^Bi, A,Aj and Bi c Bi(iJ,k,l=l,2;i<j;k<l) 
in the current adult population, respectively. Let DAiAjBtB; be 
the zygotic LD between genotypes AjAj at locus A and B^B/nt 
locus B for two-locus genotype AjAjB^Bi, i.e. DA t AjB k Bi = 
P AtAjB k B k P AiAjP B k Bj ■ There are eight zygotic LDs in total, 
but only four of them are independent since the following 

2 2 
constraints hold: ^ D Al AjB k B,= D Al A,B k B, =0. This 

ij=l;i<j k,l = \;k<l 

2 2 

leads to ^ Da^^Bj = Da 1 a 2 b,b 1 and ^ D Al A;B k B,= 

ij=l k,l=\;k<l 

2 

Yl D A,A j B Vk = -D Al A 2 B l B 1 [5\, [6]. 
ij=\\i<j 

In the continent population, let QA,AjB k Bi , Qa,Aj , and Q Bk B, be 
the frequencies of genotypes AiAjB^B/, AiAj and 
BkBi(iJ,k,l=\,2;i<j;k<l), respectively. Let Qa^J^ab, and 
D 'AiAjB k B,be the frequency of gamete^/i?^, the gametic LD, and 
the zygotic LD in the current adults, respectively. Similar 
constraints for zygotic LDs to the case in the island population 
hold as well. All zygotic and gametic LDs are assumed to be 
constant in the continent population. 

Let WA,B k (P) and WA,B k (0) be the fitness of gamete AjB^ in pollen 
and ovules, respectively. The average fitness in pollen and ovules, 
denoted by Wp and Wo, respectively, can be expressed as 
22 22 

WP = 53 53 W *Bk(P)PA ( B k and W 0=J2Y1 W ^B k (0)P* AlBk where 
i=l 1=1 1=1 1=1 

PA,B k ( = (l ~ m P)P*A,B k + m pQ*A,B K ) is tne gametic frequency after 
pollen flow. The gametic frequencies in ovules remain unaltered 
since ovules do not move after pollen flow. Let WA,AjB k B, be the 
fitness for AiAjBkB/ (ij = l,2,i<j;k,l=l,2,k<l). The average 
fitness in the sporophyte stage, w, can be calculated by 
2222 

W= Y. Y.Y. WA ' A J B ^P*A i A ] B kBl w]citrt P%A jBk B, « the ge- 

i=l j=lj<ik=\ 1=1 
notypic frequency after seed flow. Following the plant life cycle, 
the genotypic frequency after selection in the sporophyte 
stage, denoted by P% AjBkBl { = w At AjB k B,P* A . AjBtBi /w ) for A,AiB k Bi 

(ij,k,l= 1,2; i<j\k<l), is derived in Appendix SI. 

After genetic drift, the numbers of distinct genotypes follow a 
multinomial distribution. Here, the genetic sampling of jVbreeding 
individuals (an effective population size) is analogous in technique 
to but different in biological meaning from the statistical sampling 
of JV individuals [9]. Allelic and genotypic frequencies fluctuate but 
eventually reach steady-state distributions under the joint effects of 
migration, selection, and genetic drift. Gametic and zygotic LDs 
can eventually reach steady-state distributions as well. Since the 
probability density functions (pdf) of gametic and zygotic LDs are 
difficult to derive, their expectations can be indirecdy evaluated 
through multiple independent simulations. 

Genetic drift at each generation can cause the associations 
between gametic and zygotic LDs or between different zygotic 
LDs due to their sharing of some alleles or genotypes. There are 
four types of covariances between gametic and zygotic LDs, 
cov(D Al B l ,D Al AjB k B,)(ij=l,2J<j;kJ=l,2,k<l), and six types 
of covariances between distinct zygotic LDs. Note that other high- 
order LDs, such as trigenic and composite LDs [9], are not 
examined here although they can be calculated with more 
complicated analyses. Fisher's delta method is used to approximate 
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the covariances between gametic and zygotic LDs or between 
different zygotic LDs (high-order LDs) ([9], pi 18), [42]. 

For example, the covariance between Z)^,a,and Da^a^iBi or 
between D^^^and Da 1 a 2 b 1 b 2 is derived as: 

cov(Da 1 b 1 ,Da 1 a 1 b 1 b 1 ) = 

X ((! ~Pa x a x -Pb\b x )0- ~P\ -p\)p7 x a x b x b x 

-P*B* lBl ( l ~ 2 Pb\ -P*A\)PA\A l B l B 2 / 2 +P*B\PB* l B i P*A\A l B 2 B 2 0) 

-Pa\a^-Pb\-1Pa\)Pa\a 2 b 1 b 1 +Pa ] Pa 1 a 1 Pa 2 a 2 b 1 b 1 ' 

~(Pa\b x -IPA^B^iPA^B^ -ZPA^PB^)) 



COv(D A 1 A 1 B l B 1 ,E > A 1 A 2 B 1 B 2 ) — 




(P*a\a 2 b 1 b 2 2 P*a 1 a 2 Pb 1 b 2 )(Pa\a 1 b 1 b ] ^^a^b^b^J' 

where ^**is the allelic or genotypic frequency after selection in the 
sporophyte stage but before genetic drift. Other high-order 
covariances can be derived in a similar way. These covariances 
are needed for calculating the expectations of zygotic LDs 
described in the section of Analytical Theory [10]. Note that the 
above covariances produced by genetic drift are conditional on the 
allelic and genotypic frequencies />**before genetic drift, i.e. the 
expectation on the basis of genetic drift (see -E,$in Appendix S4). 
These high-order covariances can reach steady-state distributions, 
and their means, e.g.,E(cov(D Al A l B l B l ,D Al A 2 B l B 2 )), can be 
calculated in theory according to their joint probability density 
distribution (the expectation E is based on the pdf,# > described in 
the section of Analytical Theory). Similarly, multiple independent 
simulations can be used to evaluate the expectations of these high- 
order LDs. 

Note that the above general model can reduce to specific models 
with different numbers of evolutionary forces (e.g., the model with 
a single evolutionary force). Also, I only concentrate on the 
covariances between allelic frequencies, or between genotypic 
frequencies, or between gametic and zygotic LDs, or between 
different zygotic LDs. The expectations of their normalized values, 
like the square of normalized gametic LD, [8] or Lewontin's 
D'[43], are difficult to derive under genetic drift effects [44], and 
hence are not explored further. 

Monte Carlo Simulations 

MC simulations are used to examine how different evolutionary 
forces change zygotic LDs and other types of covariances in the 
plant species of a mixed mating system. Suppose that the island 
population initially has the same genetic composition as the 
continent population. For simplicity, notation for the alleles and 
subscripts in the above exact model is changed as A for^4i, a for 
A2, Bior B\, and b for B2. Simulations are conducted according to 
the plant life cycle. Given a set of parameters, including the 
genotypic frequencies in migrants and in the initial island 
population, the selection coefficients and the effective population 
size, the genotypic frequencies before genetic drift are calculated 
from Eqs. (Al) ~ (A5) in Appendix SI. For the genetic drift, 



random samples are generated using the genotypic frequencies 
that follow a multinomial distribution. Random numbers with 
uniform distribution within (0, 1) for sampling purpose are 
generated using the routine of Press et al. ([45], pp. 210-21 1). Ten 
thousand independent simulation runs are conducted for each 
case. The replicates are used to estimate means and standard 
deviations of zygotic LDs and other covariances. 

Mating System. To examine the effects of mating system, I fix 
all other parameters except the selfing rate a. Here, gametic and 
zygotic LDs are not further decomposed into the components of 
identity (inbreeding in recent ancestry) and non-identity disequi- 
libria [3], [4], and hence the interaction between selfing and 
genetic drift is unnecessarily specified. Simulations confirm that 
gametic and zygotic LDs and other covariances can reach steady- 
state distributions. Note that the parameter settings in all 
numerical examples are arbitrary as long as these parameters 
are biologically meaningful. Figure la (a coupling linkage phase, 
Dab>0) shows that the steady-state gametic and zygotic LDs have 
different patterns although they exhibit non monotonic trends with 
the selfing rate. Their standard deviations also exhibit non 
monotonic trends with a (Figure lb). Thus, gametic and zygotic 
LDs are not a linear function of a, similar to the result in a 
cytonuclear system [46]. An overlap between the steady-state 
E(D Aa BB ) and E{Daabi,) is expected when the two loci initially have 
the same settings in selection coefficients and genotypic frequen- 
cies. There are the same patterns between the steady-state 
E(cov(D A b, DAABbj) and E(cov(D A b, DaoBb)), or between the steady- 
state EicovlDjjBB, Daobb)) and E{cov(Paabb, D AaB B)), or between the 
steady-state E(cov(D AaBB , D AaBb )) and E{coiifiAABh D AaBi )) in this 
example. Selfing increases homozygosity but reduces heterozygos- 
ity, resulting in different patterns among gametic and zygotic LDs. 
The steady-state E{Dab) rnay be smaller than the steady-state 
expectations of some zygotic LDs in a predominant selfing species 
(e.g., Elfins ) in Figure la). 

The steady-state covariances between gametic and zygotic LDs 
(Figure lc) or between distinct zygotic LDs (Figure le) exhibit a 
monotonic pattern with a. Selfing facilitates the covariances 
between gametic and zygotic LDs for the genotypes with 
homozygotes at one locus (the steady-state E^oi^Djm, DAABbj) and 
E(cov(D A b, Dahbb))) or at two loci (the steady-state E(cov(D A b, 
Daabb)), but reduces the steady-state E(cov[Dab, D AaBt )). Their 
standard deviations exhibit different patterns with the selfing rate 
(Figure Id). Selfing also facilitates the covariances of zygotic LDs 
between the genotypes sharing one homozygote (the steady-state 
E(cov(Daabb, DaoBb)) and E(cov(Daabb, Daabi^I)) or sharing one 
heterozygote (the steady-state E(cov[D AaB B> D AaBt )) and E(cov(Daabi>, 
DAaBb))), but reduces the covariances of zygotic LDs between the 
genotypes without any common genotypes (the steady-state 
E^oviD^B, D AaBb )) and E(cov{D AaB B, D AABi ))). The standard 
deviations for these high-order LDs are stable with the selfing 
rate except their slight increases at the complete selfing (no effects 
from pollen flow at 0C= 1; Figure If). 

The steady-state E( D AB ) and E(D AAB b) exhibit different patterns 
with the selfing rate between the repulsion (Ats<0) and coupling 
(Ais>0 ) linkage phases although they have similar patterns in 
each linkage phase. The steady-state E(D AaB b) and E(D Aa BB) (or 
E(D A ABb)) display similar patterns withorin each linkage-phase. 
Patterns are also similar between two linkage-phase cases for the 
steady-state E{cov(Dab, AtAEffi)), but not for other three covariances 
between gametic and zygotic LDs. Selfing reduces the absolute 
steady-state E(cov(D AB , D AaB bj) in each linkage phase. All 
covariances between different zygotic LDs have the same 
responding patterns to the selfing ratecdn each linkage phase 
(data not shown here). 
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Figure 1. Effects of selfing on the steady-state gametic and zygotic LDs and other types of covariances. Average steady-state gametic 
and zygotic LDs (a) and their standard deviations (b); average steady-state covariances between gametic and zygotic LDs (c) and their standard 
deviations (d); and average steady-state covariances between distinct zygotic LDs (e) and their standard deviations (f). Results are obtained from 
10000 independent simulation runs. Parameter settings are the recombination rate = 5%, the immigration rate of pollen m P = 0.08 and seeds 
m s = 0.04, the effective population size = 50, the fitness in the gametophyte stage (pollen and ovules) w AB = 1, w„t, = 0.98, w oB = 0.98, w ab = 0.96, and 
the fitness in the sporophyte stage w AABB = 1 , w AABb = w AaBB = 0.98, w Mbb = w AaBb = w aaBB = 0.96, m Aabb = w aaBb = 0.94, and w aabb = 0.92. The genotypic 
frequencies in the continent and initial island populations are 0.1225 for AABB, 0.105 for AABb, 0.0225 for A Abb, 0.105 for AaBB, 0.245 for AB/ab, 0.045 
for ab/aB, 0.105 for Aabb, 0.0225 for aaBB, 0.105 for aaBb, and 0.1225 for aabb. 
doi:1 0.1 371 /journal. pone.0080538.g001 
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The above examples indicate that plants with distinct mating 
systems have different zygotic LDs and other covariances in a local 
population even under the same impacts of immigration. Both 
zygotic and gametic LDs are sensitive to the pattern of mating 
system. Predominant outcrossing species have weaker covariances 
between zygotic and gametic LDs, but stronger covariances 
between distinct zygotic LDs than do the predominant selfing 
species. 

Seed and Pollen Flow. To examine the effects of pollen (or 
seed) flow, I fix all other parameters except the migration rate of 
pollen (or seeds). Figure 2a shows that the steady-state £(Aib) 
slightly increases withmp in a coupling linkage phase (D^g>0). 
The steady-state Elfljm^ ) and E(D AllBb ) slightly decrease with mp, 
while the steady-state E(D AaBB ) and E(D AABb ) (negative) slightly 
increase with mj>(Figure 2a). The steady-state ElcovlDj^, D AABb )), 
E(cov(D AB , D AaBB )), and E(cov(D AB> D AaBb )) slightly increase with mp 
while the steady-state E(cov(D AB , A^is^)) decreases with mp(Figure 
2c). All covariances between different zygotic LDs slightly decrease 
withm/>(Figure 2e). All standard deviations slighdy decrease with 
»ip(Figures 2b, d, and f). 

Seed flow has greater effects than pollen flow on zygotic LDs 
and other covariances (Figure 3; the same parameter settings as in 
Figure 2 except the different migration rates of seeds and pollen). 
The steady-state £(Aig) changes faster withm^than any steady- 
state zygotic LDs. Generally, the steady-state zygotic LDs and the 
covariances between different zygotic LDs or between gametic and 
zygotic LDs do not monotonically change with nig. The steady- 
state E( D^jj^ and E(D AaBb ) slightly increase asOT^approaches the 
value of selection coefficient, and then slightly decrease afterwards 
(Figure 3a). Similar patterns exist for the change of the steady-state 
£(«)i(AiB, Djubb)), E^Dju;, D AaBb )), E(cov(D AABB , D AaBb )), and 
E(cov(D AaBB , Djytgh)) with /^(Figures 3c and e). To the contrary, 
the steady-state E(D AaBB ) and E(D AABi ) slightly decrease as 
Wjapproaches the value of selection coefficient, and then slightly 
increase afterwards. Similar patterns exist for the change of the 
steady-state E(cov(D AB , D AaBB )), £[co»(Aib, Amj»))> E{cov{D aabb , 
E>AoBb)), Eicov^D^jj, DjUBt)), E[cotiD AaBB , D AaBb )), and E^oriD^,,, 
E> Aa Bb)) (Figures 3b and c). The same pattern exists for the 
covariances between zygotic LDs for the genotypes with a 
common genotype at one locus, or for the genotypes without 
any common genotype at each locus. All standard deviations 
gradually decrease with ms (Figures 3b, d, and f). 

These examples indicate that a local plant population generally 
exhibits robust responses to the impacts of immigration of pollen 
and seeds in terms of zygotic LDs, or the covariances between 
gametic and zygotic LDs, or the covariances between distinct 
zygotic LDs. Seed and pollen flow have small effects on high-order 
LDs in a local population. 

Genetic Drift. To examine the effects of genetic drift, I fix all 
other parameters except the effective population size (jV). Figure 4 
shows the results for a predominant outcrossing species (a = 5%). 
The steady-state E(D AB ) and E{D AABB ) slightly increase with JV. 
The steady-state E(D AaBB ), E(D AABb ), and E(D AaBb ) (genotypes with 
heterozygote at one locus or two loci) slighdy decrease as the 
effective population size increases (Figure 4a). The steady-state 
E(cov(D AB> Ambb)) and E(cov(D AB , D AaBb )) gradually reduce to zero 
as JV increases. To the contrary, the steady-state E{cov(D AB , A^iz;*)) 
and E(cov(D AB , JJ AaBB j) gradually increase to zero as JV increases 
(Figure 4c). The steady-state E^oiiD^s, D AaBh )) and E[cov(D AaBB , 
E>aab£)) gradually reduce to zero with JV, while other steady-state 
E^oviDAASB, D AaBB )), E( covings, D AABb )) ! E(cov(D AaBB , D AaBb )), and 
E(cov{D AABb , D AaBb j) gradually increase to zero with JV (Figure 4e). It 
is clear that covariances between gametic and zygotic LDs or 
between different zygotic LDs are more sensitive than gametic and 



zygotic LDs to the genetic drift effects. All standard deviations 
gradually decrease with JV (Figures 4b, d, and f). 

The examples indicate that a small local population can affect 
zygotic LDs, and has large effects on the covariances between 
gametic and zygotic LDs or between distinct zygotic LDs. These 
high-order covariances are more informative than gametic LD in 
signaling the effects of population demographic dynamics. 

Selection. To assess the effects of linear-additive selection, I 
examine three selection schemes: gametic selection only, zygotic 
selection only, and both gametic and zygotic selection. Table 1 
shows a comparison in the steady-state zygotic and gametic LDs 
and other types of covariances. The steady-state -E^D^g) slighdy 
decreases while the absolute steady-state zygotic LDs and other 
types of covariances increase from the case of gametic selection 
only to the case of zygotic selection only, and to the case of joint 
selection. The examples indicate that cumulative selection can 
enhance zygotic LDs and other covariances in the linear additive- 
viability model (Table 1). 

To assess the effects of epistatic selection, I use Dobzhansky- 
Muller's incompatibility model [27], [28], [47] as an example to 
demonstrate how epistatic selection in the sporophyte stage affects 
gametic and zygotic LDs. Three cases with different extents of 
epistatic selection are examined. Selection in the gametophyte 
stage is excluded in each case. In Case I, the genotypic fitness is set 
as w AABB = w aatb = 1, w AaBB = w, mBb = 0.99, w AABb = w Aabb -0.99, 
w aaBB =0.98, 1X^444 = 0.98, and w AaBb = 0.98. In Case II, the 
genotypic fitness is set as w AABB =w aabb = 1, w AaBB = w aaBb = 0.99, 
w AABb = w Aabb = 0.5, w aaBB = 0.98, w AAbi = 0.5, and w AllBb = Q.5. In 
Case III, the genotypic fitness is set as w AABB = u> aabb = 1 , 
«>AaBB = WaaBb = 0.99, wjum = w Aaib = 0.1, w aaBB = 0.98, Wjuu = 0.1, 
and w AaBb = 0.l. These three cases are the same as matrices (13), 
(14), and (15) of Gavrilets [48], respectively. In these settings, 
alleles A and b have a progressively negative interaction on fitness 
(incompatible background interactions) from Cases I to III. 

Results indicate that epistatic selection can change the relative 
gametic and zygotic LDs (Table 2). The steady-state frequency of 
allele B increases while the steady-state frequency of allele A 
decreases from Cases I to III. The steady-state E(D AaBB ) and 
absolute steady-state E(D AaBb ) become greater than the steady-state 
E(D AB ) in Case III. The steady-state E(D AB ), E(D AABB ), and 
E(D AaBb ) decrease while the steady-state E(D AaBB ) and E(D AABb ) 
increase from Cases I to III. Epistatic selection also changes the 
covariances between gametic and zygotic LDs or between different 
zygotic LDs. The steady-state E(cov(D AB , D^n^), E(cov(D AB> 
E> AABb j), and E(cov{D AB , D AaBb )) decrease while the steady-state 
E(cov(D AB> D AaBB j) increases from Cases I to III. The steady-state 
EicoviD^j;, Aw;*)), E[cov{D AABB , D AaBb )), and E{cov(D AaBB , AlW) 
decrease while the steady-state E(cov{D AABB , D AaBB )), E(cov(D AaBB , 
E> Aa Bb)), and E(cov(D AABb , D AaBb )) increase from Cases I to III. 

The above examples indicate that zygotic and gametic LDs 
have different responding patterns to natural selection. The 
cumulative selection can enhance zygotic LDs and other 
covariances in the additive-viability selection model. One striking 
result is that epistatic selection at the diploid level can produce 
zygotic LDs that are greater than or comparable to gametic LD. 
This pattern can be used to detect the epistatic selection process in 
natural populations. 

Analytical Theory 

To further understand the evolution of zygotic LDs, I derive the 
analytical theory in a linear-additive-viability model with weak 
selection and random mating (a = 0). The gametic fitness in pollen 
and ovules is decomposed as w AiBk (p) = 1 +s Ai (p) + s Bt (P) and 
WA l B k (0) = \+s Ai (o)+SB k (0) where s A ,(p) and s Ai(0) are the 
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Figure 2. Effects of pollen flow on the steady-state gametic and zygotic LDs and other types of covariances. Average steady-state 
gametic and zygotic LDs (a) and their standard deviations (b); average steady-state covariances between gametic and zygotic LDs (c) and their 
standard deviations (d); and average steady-state covariances between distinct zygotic LDs (e) and their standard deviations (f). Results are obtained 
from 10000 independent simulation runs. Parameter settings are the selfing rate =5%, the recombination rate = 5%, the effective population 
size = 50, the immigration rate of seeds m s = 0.04, and the fitness in the gametophyte stage (pollen and ovules) w AB ="\, w Ab = 0.98, iv„ e = 0.98, 
w ab = 0.96, and the fitness in the sporophyte stage w AABB = 1 , w MBb = w AaBB = 0.98, w AAbb = w AaBb = w aaBB = 0.96, w Aabb = w aaBb = 0.94, and w aabb = 0.92. 
The genotypic frequencies in the continent and initial island populations are 0.1225 for AABB, 0.105 for AABb, 0.0225 for AAbb, 0.105 for AaBB, 0.245 
for AB/ab, 0.045 for ab/aB, 0.105 for Aabb, 0.0225 for aaBB, 0.105 for aaBb, and 0.1225 for aabb. 
doi:1 0.1 371 /journal.pone.0080538.g002 



selection coefficients for allele A, in pollen and ovules, respectively; 
s B t (P) an d SB k (0) ar e the selection coefficients for allele in pollen 
and ovules, respectively. The genotypic fitness in the sporophyte 



stage is expressed as WAiAjB k B, = 1 +SAtAj +SB k B t where ^,-^.and 
.y^ t £,are the selection coefficients for genotypes AjAj and B^B/, 
respectively. 
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Figure 3. Effects of seed flow on the steady-state gametic and zygotic LDs and other types of covariances. Average steady-state 
gametic and zygotic LDs (a) and their standard deviations (b); average steady-state covariances between gametic and zygotic LDs (c) and their 
standard deviations (d); and average steady-state covariances between distinct zygotic LDs (e) and their standard deviations (f). Results are obtained 
from 10000 independent simulation runs. Parameter settings are the selfing rate = 5%, the effective population size = 50, the immigration rate of 
pollen m P = 0.04, and the fitness in the gametophyte stage (pollen and ovules) w AB ="\, w Ab = 0.9S, w og = 0.98, w ab = 0.96, and the fitness in the 
sporophyte stage w MBB = 1 , w MBb = w AaBB = 0.98, w AAbb = w AaBb = w aaBB = 0.96, w Aabb = w aaBb = 0.94, and w aabb = 0.92. The genotypic frequencies in the 
continent and initial island populations are 0.1225 forces, 0.105 for AABb, 0.0225 for AAbb, 0.105 for/loee, 0.245 for AB/ab, 0.045 for ab/aB, 0.105 for 
Aabb, 0.0225 for aaBB, 0.105 for aaBb, and 0.1225 for aabb. 
doi:1 0.1 371 /journal.pone.0080538.g003 

With the weak selection, all items containing the second or items containing the second or higher orders of the migration rate 
higher orders of selection coefficients are neglected. The immi- m|, or m s m P , or higher orders) or the products of the 

gration rates of seeds and pollen are assumed to be small. The migration rate with selection coefficients (sm P or sm s ) are 
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Figure 4. Genetic drift effects on the steady-state gametic and zygotic LDs and other types of covariances. Average steady-state 
gametic and zygotic LDs (a) and their standard deviations (b); average steady-state covariances between gametic and zygotic LDs (c) and their 
standard deviations (d); and average steady-state covariances between distinct zygotic LDs (e) and their standard deviations (f). Results are obtained 
from 10000 independent simulation runs. Parameter settings are the selfing rate = 5%, the immigration rate of seeds m s = 0.04 and pollen m P = 0.08, 
and the fitness in the gametophyte stage (pollen and ovules) w AB =\, 1^ = 0.98, w oB = 0.98, w ab = 0.96, and the fitness in the sporophyte stage 
Waabb=1, w MBb = w AaBB = 0.98, w AAbb = w AaBb = w aaBB = 0.96, w Aabb = w aaBb = 0.94, and w aabb = 0.92. The genotypic frequencies in the continent and 
initial island populations are 0.1225 for AABB, 0.105 for A4B6, 0.0225 for AAbb, 0.105 for AaBB, 0.245 for AB/ab, 0.045 for ab/aB, 0.105 for Aabb, 0.0225 
for aaBB, 0.105 for aaBb, and 0.1225 for aabb. 
doi:1 0.1 371 /journal.pone.0080538.g004 



neglected. Again, notation for the alleles and subscripts in the 
exact model is changed as A for A\, a for A%, B for B\, and b for 
Bi. Selection coefficients are set as SAi(p) = Sa 1 (O) = 0> 



$B,(P) — SBjiO) =0, SA 2 {P)——S a p, S Al ( 0 )——S a0 , S B2 (p)— —S h p, 

and Sb 2 (0) = —Sho- Alleles a and b are maladaptive in the island 
population. Let h a and ht are the degrees of dominance at loci A 



PLOS ONE | www.plosone.org 



8 



November 2013 | Volume 8 | Issue 11 | e80538 



Evolution of Zygotic LD 



Table 1. Effects of selection in the gametophyte stage, the sporophyte stage, and in both stages on the steady-state gametic and 
zygotic LDs and other high-order covariances* 





Gametic selection 


Zygotic selection 


Gametic and zygotic selection 


Pa 


U.J/OU— U. I ZZ I 


n ^rch + n 1917 

U.J03J— U. 1 Z. 1 / 


n ft^9i +n 1 1 ^ 

U.OjZ 1 _U. 1 1 JJ 


PB 


n E iR99 + n 1 99zl 

U.JOZZ _ U. 1 ZZH 


fl ^01 Q + f) 1 9flQ 


n ft^^Q+n 1131 

U.UJJ3 — U. I I J I 


U AB 


U.U / JJ — U.UtJ/ 


U.U / H-U _ U.utJJ 


U.U/ jU_U.UtDJ 


&AABB 


0 061 9 + 0 0AQ5 


0 069Q + 0 0A9Q 


0 0799+f) 0519 

U.U/ Z.Z. — U.UJ 1 _7 


u AaBB 


n rn^R+n rtA97 

— U.UjjO — U.UtZ / 


U.U J J J _ U.UtJJ 


U.UH-DO _ U.UH- JO 


&AaABb 


fl 0^49 + 0 0A9R 

U.U JtZ —U.UtZO 


U.UJJO — U.UtJ 1 


U.Ut/ I —U.UtJO 


U Aa Bb 


0 fiA37 + n OAftR 


U.UtJO _ U.UH-O/ 


fl flAQ9 +n fiA71 


cov(D AB ,D AABB ) 


i fiiQvin~ 4 +i mivin -4 

J.DI7AIU _ 1. IU 1 A IU 


^ fi7Qv1fl~ 4 +1 1f17 y1D~ 4 

J.D/7 a IU _ I . IU/ A I U 


An^ c ivin _4 +1 riQQvin -4 

H.UJJ A 1 U _ 1 .U7J A 1 U 


CO\(D AB ,D AaBB ) 


1 Ri3vin _4 +1 ^ai vin -4 

I.OI J A IU _ I .0*+ 1 A 1 U 


1 RR7vin _4 +i ^Wvin -4 

I .00/ A IU _ I.jjUa IU 


9 3ARvin _4 +1 39Qy1fl -4 

Z.3H-0 A I U _ I .JZ7 A 1 U 


cov(D AB ,D AABh ) 


-1.840x10" 4 ±1.329x10" 4 


-1.904x10" 4 ±1.338x10" 4 


-2.368x1 (T 4 ±1.31 3 xlCT 4 


cov(D AB ,D AaBh ) 


2.691 x10" 4 ±1.615x10" 4 


2.689x10" 4 ±1.616x10" 4 


2.880x10~ 4 ±1.508x10~ 4 


cov(D AABB ,D AaBB ) 


-6.525 x10" 4 ±2.086x10" 4 


-6.653 x 1 0" 4 ±2.064 x 1 0" 4 


-7.266 x1(T 4 ±1. 877 x10~ 4 


cov(D AABB ,D AABh ) 


-6.550 x10" 4 ±2.088x10" 4 


-6.665 x 1 0" 4 ±2.056 x 1 0" 4 


-7.280 x10~ 4 ±1. 874 x10~ 4 


cov(D AABB ,D AaBb ) 


5.107x10" 4 ±2.043x10" 4 


5.256x10" 4 ±2.025x10" 4 


6.088x10~ 4 ±1.870x10~ 4 


cov(D AaBB ,D AABb ) 


5.500x10" 4 ±2.083x10" 4 


5.636x10" 4 ±2.059x10" 4 


6.428x10~ 4 ±1.855x10~ 4 


cov(D AllBB ,D AaBh ) 


-6.816x10" 4 ±2.344x10" 4 


-6.945 x10" 4 ±2.297x10" 4 


-7.631 x1(T 4 ±2.035x1(r 4 


co\(D AABh ,D AtlBf} ) 


-6.723x10" 4 ±2.338x10" 4 


-6.892 x 1 0" 4 ±2.31 6 x 1 0" 4 


-7.562x10~ 4 ±2.051 x10~ 4 



Three selection schemes are: w AB = 1 ,w Ab = 0.98, w a8 = 0.98, and w ab = 0.96 for gametic selection only; w AABB ='\, w AABb - w acbb - 0-98, w^ bb = w/^ aeb = w aa8B = 0.96, 
W/iott = w m8b = 0-94, and w aabb = 0.92 for zygotic selection only; and w AB = 1 ,w Ab = 0.98, w oB = 0.98, w ab = 0.96, iv^,, = 1 , = w AaBg = 0.98, w M66 = w AaBb = w aaBB = 0.96, 
w Aabb = w aoefe = 0-94, and = 0.92 for both gametic and zygotic selection. Other parameter settings are the recombination rate = 5%, the immigration rate of pollen 
m P = 0.08 and seeds m s = 0.04, and the effective population size = 50. The genotypic frequencies in the continent and initial island populations are 0.1225 for AABB, 0.105 
for AABb, 0.0225 for AAbb, 0.1 05 for AaBB, 0.245 for AB/ab, 0.045 for ab/aB, 0.1 05 for Aabb, 0.0225 for aaBB, 0.1 05 for aaBb, and 0.1 225 for aabb. The steady-state results 
(mean ± S d ) are obtained from 10000 independent simulation runs. 
doi:1 0.1 371 /joumal.pone.0080538.t001 



and B, respectively. Selection coefficients for genotypes are set as 
•5/41^1=0, s A x A 1 = —h„s a , and SA 1 A 2 = ~ s a for locus A; and 
SBtBi =0, Sb { b 2 = —hbSb, and sb 2 b 2 = — s b f° r locus B. 

From Eqs. (Al) ~ (A5) in Appendix SI, the deterministic 
changes in allelic frequency (ApA andAps), gametic (ADab) and 
four independent zygotic LDs (ADaabb, AD^ABb, ADauBB, and 
ADAuBb), can be derived. Other functions of zygotic LDs can be 
calculated once the four independent zygotic LDs are available. 
After genetic drift, the means for the per-generation changes in 
allelic frequency, gametic and zygotic LDs, can be derived using 
the conventional approach [49] (Appendix S2). Note that one 
additional factor ( 1 — f)timesZ>^_g in the formulae in Appendix S2 
is because Dab is termed from the preceding adults in a plant life 
cycle (one generation difference between adults and pollen and 
ovules; [41]). 

Let <P(pA,PB,D A B,DAABB,DAABb,DAaBB,DAaBb)be the steady- 
state pdf at the two linked loci so that <Pdp AdpsdD ab^D aabbA 
D AABbdD AaBBdD AaBb represents the expected number of two loci 
having the allele frequencies, gametic and zygotic LDs within the 
intervals (p A ,PA+dp A ), (pB,PB + dp B ), ■ and (D Aa Bb, D Aa Bb + 
dDA,,Bb ), respectively. Expectation of each individual variable can 
be calculated in theory from pdf 0. For instance, an expectation of 

gametic LD can be obtained byE(DAB)= ■■■ QDABdpAdpsd 

D aabb- -dD AaBb- For a stationary distribution of a function of 
seven variables g(pA,PB,D A B,DAABB,DAABb,D A aBB,DA a Bb), the 
Kolmogorov backward equation can be derived in the following 
expression [12], [50]: 



Q = E{ M(d p .)-^+M(S p „) 



('Pa 



SpB 



8g 



v AB 'dD 



+cov(d PA ,S PB ) 



dg 



dp A 8p B 



A " Bt ' 0 D A aBb 

+ cow(6.,,,3 D , R ) 



c 1 



8 2 g 



SpaSDab 



(3) 



■+^AaB B MaBb) 8DAaBBdDAaBb 



2 " A dpi 2 " B dp\ 2 AaBb dD\ aBh 



Notation E in Eq. (3) means expectation with respect to pdf 0, 
the same meaning as in the preceding section except that its 
calculation is based on numerical simulations. 

In Eq. (3), there are seven items with the average change 
coefficientsM((5 ), seven items with the variance coefficients V(8), 
and twenty-one items with the covariance coefficients. Appendix 
S3 gives the expressions for the variances of per-generation 
changes in allelic frequency, gametic and zygotic LDs, and all 
possible covariances among these per-generation changes. 

With the diffusion model, the expectations of zygotic LDs and 
the covariances between gametic and zygotic LDs or between 
different zygotic LDs can be calculated in theory. However, the 
algebraic deduction remains complicated when the joint effects of 
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Table 2. Effects of Dobzhansky-type epistatic selection 
covariances* 


on the steady-state gametic 


and zygotic LDs and other high-order 






Case 1 


Case II 


Case III 


PA 


0.5068± 0.1280 


0.3090±0. 2127 


0.2673±0.2343 


PB 


0.5078±0.1292 


0.7573±0.1937 


0.8170±0.2009 


Dab 


0.0788±0.0510 


0.0329±0.0218 


0.0165±0.0113 


Daabb 


0.0533 ±0.0462 


0.0219±0.0236 


0.0100±0.0112 


DauBB 


-0.0235±0.0394 


0.0257±0.0337 


0.0310±0.0248 


DAaABb 


-0.0231 ±0.0397 


-0.0172±0.0205 


-0.0079±0.0095 


DAaBb 


0.0436±0.0474 


-0.0 120 ±0.03 12 


-0.0245 ±0.02 15 


cov(D A b,Daabb) 


3.217x10~ 4 ±1.149x10~ 4 


1.368x10~ 4 ±1.195x10~ 4 


0.556x10~ 4 ±7.466x10~ 4 


cov(D AB ,D AllBB ) 


-1.384x10~ 4 ±1.353x10~ 4 


0.6650 x 10~ 4 ± 1.1 58 x10~ 4 


8.250x10~ 4 ±6.508x10~ 4 


COV(D AB ,D A ABb) 


-1 . 399 x10~ 4 ± 1.341 x10~ 4 


-0.9688 x 1 0~ 4 ±0.8927 x 1 0~ 4 


-0.377 x 1 0~ 4 ±0.467 x 1 0~ 4 


co\(D A B,D AaB b) 


2.692 x10~ 4 ±1. 642 x10~ 4 


0.1426x10~ 4 ±1.194x10~ 4 


-0.499 x 1 0~ 4 ±0.505 x 1 0~ 4 


00\(D AABB^D AuBB) 


-5.511 x10~ 4 ±2.298x10~ 4 


-0.928 x10~ 4 ±1. 275 x10~ 4 


0.1252x10~ 4 ±0.537x10~ 4 


COv(D AA BB,DAABb) 


-5.518x10~ 4 ±2.289x10~ 4 


-1.987x10~ 4 ±1.753x10~ 4 


-0.6276 x 1 0~ 4 ±0.823 x 1 0~ 4 


COV(D AABB ,D AaBb ) 


3.942 x10~ 4 ±2.1 05 x10~ 4 


0.8561 x10~ 4 ±1.252x10~ 4 


-0.092 x 1 0~ 4 ±0.486 x 1 0~ 4 


cov(D AaBB ,D AABh ) 


4.386x10~ 4 ±2.212x10~ 4 


0.8903 x10~ 4 ±1. 284 x10~ 4 


-0.091 x10~ 4 ±0.489x10~ 4 


COv(D AaBB ,D AaBb ) 


-5.628 x 1 0~ 4 ±2.536 x 1 0~ 4 


-5.495 x 1 0~ 4 ±2.362 x 1 0~ 4 


-2.307 x10~ 4 ±1. 544 x10~ 4 


cov(D AABb ,D At i Bh ) 


-5.611 x10~ 4 ±2.513x10~ 4 


-0.7594x10~ 4 ±1.251 x10~ 4 


0.128x10~ 4 ±0.4867x10~ 4 



•Three selection schemes are: w MBB = w aabb = 1 , w AaBB = w aaBb = 0.99, w MBb = w Aabb = 0.99, w aaBB = 0.98, w Mbb = 0.98, and w AaBb = 0.98 for Case I; w MBB = w aabb = 1 , 
Waobb = WaoBt, = 0.99, w MBb = w Aabb = 0.5, w oo8B = 0.98, w Mbb = 0.5, and ^^ = 0.5 for Case II; w MBB = w aabb = 1 , w AaBB = w aaBb = 0.99, w MBb = w Aabb = 0.1, w aaBB = 0.98, 
w AAbb = 0.'\, and w Aagfa = 0.1for Case III. Other parameter settings are the recombination rate = 5%, the immigration rate of pollen m P = 0.08 and seeds m 5 = 0.04, and the 
effective population size = 50. The genotypic frequencies in the continent and initial island populations are 0.1225 for AABB, 0.105 for AABb, 0.0225 for AAbb, 0.105 for 
AaBB, 0.245 for AB/ab, 0.045 for ab/aB, 0.105 for Aabb, 0.0225 for aaBB, 0.105 for aaBb, and 0.1225 for aabb. The steady-state results (mean ± S d ) are obtained from 
10000 independent simulation runs. 
doi:1 0.1 371 /joumal.pone.0080538.t002 

selection, migration, and genetic drift are considered. Here, I 
consider two specific cases. One case is that locus A is selective 
while locus B is neutral, with an emphasis on genetic hitchhiking 
effects [51], [52]. The other case is that both loci are neutral, with 
an emphasis on the effects of linkage distance. 

Genetic Hitchhiking. How genetic hitchhiking effects evolve 
zygotic LDs is an important issue for studying the pattern of 
genotypic diversity along chromosomes. This may provide a 
genetic basis for forming zygotic LD blocks, analogous to the 
gametic LD blocks along chromosomes. Suppose that locus A is 
mainly subject to the balance between the effects of selection and 
immigration. The genetic drift effects are negligible for locus A. 
Locus B is subject to the balance among the effects of immigration, 
genetic drift, and recombination with locus A. This consideration 
is similar to the previous studies in examining associative 
overdominance or genetic hitchhiking effects on spreading neutral 
nuclear/organelle genes [44], [53], [54], [55]. All items withs/,p, 
S/,0, and S/, are eliminated for the average per-generation changes 
in allelic frequency, gametic and zygotic LDs in the formulae in 
Appendix S2. The variances for the per-generation changes in 
allelic frequency Pa and all covariances between and gametic 
LD or between pa and zygotic LDs are removed in the formulae 
in Appendix S3, but the remaining expressions hold except that 
the steady-state Pa is known. Similarly, the items containing 
5p A ,V(5 PA ), and the covariances between 5 PA and gametic LD or 
between 5 pA and different zygotic LDs in Eq. (3) are removed. 

The steady-state equation for allelic frequency at locus A can be 
obtained by setting ApA = 0, the same as setting M(S pA ) = 0, and 
s bP = s bo= s h = Q in Appendix S2, i.e. 



0 = (m s + m P /2)(Q A -p A ) 

(4) 

+PAPa (OaP + S a0 )/2 + S a (p a ~ K{ 1 - 2p A ))) . 



The steady-state allelic frequency can be numerically calculated 
from the above cubic equation, given the condition of 0 <Pa < 1 
and other parameters. Like Ohta and Kimura [44], denote ^^or 
p a as the known frequencies calculated from Eq. (4). It can be seen 
that selection in the gametophyte and sporophyte stages is 
compounded in the case of h a - 1/2. 

To calculate the expectations of the steady-state zygotic LDs 
and other types of covariances from Appendices S2 and S3, the 
following fourteen expectations are required: E(p A B ^'){i= 0, 1, 2, .3), 
E(pl- J D AJi )(j= 0,1,2,3), E(p 2 B - k D 2 AB )(k= 0,1,2), Etff 1 ti^) 
(/=0,1), and E(D AB ). Expectations of a few low-order functions 
can be analytically derived. For instance, letting g=Ps an d g = Dab 
separately in Eq. (3), I can obtain: 



E(D )= m(\-r)D AB f5 , 

AB > r + { i- r){ f h+ i/2N-(l-p A -Q A )Ay y ' 



m , n , (\-r) 2 AD AB 

E(PB) = QB+ r + V-r)(m+l/2N-(l-p A -Q A )Ay (6) 

where m = ms + mp/2, the joint migration rate, and 
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A = (SaP + Sao)/2 + s a (p a —h a (\—2p A )), the selection component 
at locus A. Eq. (6) indicates the dependence of the allelic frequency 
at locus B on the allelic frequency at locus A. 

Substitution of g in Eq. (3) by three functions, p B , PbDab, an d 
D AB> can yield three equations for calculating E(p 2 B ), E(pbDab), 
and E(D 2 AB ): 



E( V(D AaBh )) = ^ {4p A pa(l - 2pAPa)E(( 1 - 2p B p h )p BP b) 

+2(1-2^X1 -Ap A p a ){\ -r)E{{\ -2p B ){\ -4p B p b )D AB ) (9) 
_ 16(1 -,) 3 (1 -2p A )E({\ -2p B )D\ B )- 16(1 -r) 4 E(D 4 AB )) 



an 


an 


0 


an 


a 2 2 


A23 


a 3 i 




fl 33 



'E(p 2 B ) 

E(p B D AB ) 
,E(D AB ) 




(7) 



where an = 2m + 1/2jV, ai 2 = — 2(1 — r)A,a 2 i = ^AQa — Pa), 
and g 3 = (2mgj^ - g^) -2w(l-r)A^- (1-2^4)0-'') 
/2N)E(D AB ) -pAp a E(p B )/2N. 

Expectations of the remaining nine functions can be numeri- 
cally calculated using Mathematica tool [56] by substituting g in Eq. 
(3) withpl,p B ,p 2 B D AB ,plD AB ,p B D 2 AB ,p 2 B D AB , D AB ,p B D AB , and 
D AB> respectively. These calculations are not shown here. 

With the availability of the above fourteen expectations, the 
expectations of some lower or the same order functions can be 
indirectly calculated. For instance, I can obtain 
e (Pab) =PaE(pb) + (1 - r)E(D AB ) and E(p* AB p* ah ) =pApaE(p B Pb) 
+ (1 - r)(p A E(p B D AB )+p a E(p h D AB )) + (1 - r) 2 E(D 2 AB ). The ex- 
pectation of any steady-state zygotic LD, E(D,,„), can be 
calculated by substituting g in Eq. (3) with D (one of the four 
independent zygotic LDs), resulting in E(M(Sp> ) = 0. For 
instance, E(D AaB i,)can be calculated byE(M(SD AiM )) = 0 from 
Appendix S2, i.e. 

E(D Aa Bb) = m P ( Q* AB E(p* ah ) + QabE(p* AB ) 
+ Q Ah E(p* aB ) + Qa B E(p Ah ) 

~ 2pAPa{QbE{p B ) + QBE{pb)) 
+ m S (QAaBb-2p A p a Q Bh ) 

- 2(m s + mp)(E(p AB p* ab ) + E(p Ab p* aB )) 
-2{\-r)E{D AB )/N (8) 
+ (1—0(2(1 -2p A ) + (s a p+s a0 )(\ -6p A p a ) 
+ 2s a (p 2 a (l-4p a ) 

-K(\ -Ap A Pa){\-2pA))E((\-2p B )D AB ) 
+ 2{\-r) 2 (2- s aP - s a0 + 2s a (p a {\ - 2p a ) 
-h a {\-Ap A p a ))E(D 2 AB ) 



Eqs. (5) and (8) indicate that effects of seed and pollen flow are 
compounded in generating gametic LD, but can be separated in 
generating zygotic LDs. 

The expectations of the steady-state variances of any zygotic 
LDs can be calculated using Fisher's delta method by omitting all 
items containing m /N, sJN, and higher orders. It is shown that 
these expectations can be calculated from the expectations of the 
variances of the per-generation change in zygotic LD in Appendix 
S3, i.e. E(V(D, ^)) = E(V(5d ))with a sufficient accuracy (Ap- 
pendix S4; [12]). For instance, E(V(D Aa Bb)) can be calculated 
from E(V{& Dam )), i.e. 



The expectation of any steady-state covariance between 
gametic and different zygotic LDs can be calculated using Fisher's 
delta method by omitting all items with mjN, s./N, and higher 
orders. It is also shown that this expectation can be calculated from 
the expectation of the covariance in its per-generation change in 
Appendix S3, i.e. E(cov(D AB ,D )) = E(cov(3 Dab S d )) with a 
sufficient accuracy (Appendix S4). For instance, E(cov(D AB ,D AaB b)) 
can be calculated from E(cov(Sd ab ,^D M1) ))in Appendix S3, i.e. 

E{cow(D AB ,D AaBb )) = 



(10) 



^{pAPa(\ -2p A )E(p B p h (\-2p B )) 

+ (1 - 0(1 - 2p A ) 2 E{(\ - 2 PB ) 2 D AB ) 
+ 4(1 - r)p A p a E(p B p h D AB ) 
-(\-2p A )(\-r) 2 E{(\-2 PB )D 2 AB ) 
— 4(1 — E{D AB )) 



Similarly, expectations of other covariances in Appendix S3, 
such as the covariances between different zygotic LDs, can be 
calculated in the way similar to the above deductions. Expecta- 
tions of high-order LDs, such as E(D 2 AaBh ), can be numerically 
calculated using multiple equations derived by substituting g in Eq. 
(3) with D 2 AaBh and other high-order LDs. This needs more 
extensive algebraic analyses, and is not explored further. 

Simulations confirm that the above analytical model performs 
well. For instance, consider the same parameter settings as in 
Figure 1 for the genotypic frequencies in the continent and island 
populations, a =0, m P = 0.08 and m s = 0.04, jV= 100, 
s a(0) -s a (P) — s a /2 — 0.04, and h a =0.5. Gametic and zygotic 
LDs and other covariances can reach steady-state distributions 
(~50 th generation; data not shown here), reflecting the equilibrium 
among the effects of migration, genetic drift, and genetic 
hitchhiking. All analytical predictions are distributed within the 
ranges of one-standard deviations of the simulation results (Table 
3). 

Figure 5 shows that genetic hitchhiking effects can produce 
different patterns among gametic and zygotic LDs and other 
covariances. The expected neutral allelic frequency, E(jj£), 
increases as the frequency of favorite allele A increases with the 
selection coefficient. E(D AB ) gradually decreases while ^(AwBi)) 
gradually increases with the selection coefficient (Figure 5a). 
E{D AaB b) slighdy increases while both E(D AaBB ) and £(Alib«) 
decrease with the selection coefficient. The covariances between 
gametic and zygotic LDs for the genotypes with heterozygotes at 
one locus or two loci gradually decrease with the selection 
coefficient, except E( cov{D AB , Daabb)) showing a different pattern 
(Figure 5b). The covariances between the zygotic LDs for the 
genotypes with one common genotype at a locus decrease with the 
selection coefficient, while the covariances between the zygotic 
LDs for the genotypes without a common genotype at one locus 
(EXcotiD^B, D AaBh )) and E(cov(D AaBB , D^tfj) increase with the 
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Table 3. Comparison between the simulation results and analytical model predictions under genetic hitchhiking effects* 





MC simulations 


Analytical model 


PA 


0.7165±0.0716 


0.7071 


PB 


0.5916±0.0915 


0.5509 


D AB 


0.0777±0.0355 


0.0536 


Daabb 


0.0732 ±0.0405 


0.0465 


DauBB 


-0.0533 ±0.0347 


-0.0317 


DAaABb 


-0.0367 ±0.03 17 


-0.0169 


DAaBb 


0.041 9±0.0375 


0.0202 


co\(Dab,D AA bb) 


2.011 x10~ 4 ±4.323x10~ 5 


2.01 8 x10~ 4 


cov(Dab,D Aii bb) 


-1.463x10~ 4 ±3.994x10~ 5 


-1.389x10~ 4 


COV(D A B^AABb) 


— i .UDy x i u _ j . / y j x i u 


— / . 1 oo x 1 u 


CO\(D A B,D AaB b) 


1.541 x10~ 4 ±6.185x10~ 5 


1.032x10~ 4 


CO\(D A ABB,D A aBB) 


-3.918x10~ 4 ±7.381 x10~ 5 


-4.125x10~ 4 


cov(D A abb,£>aabi,) 


-3.496x10~ 4 ±8.167x10~ 5 


-3.514x10~ 4 


CO\(D A ABB,D Aa Bb) 


3.140x10~ 4 ±7.641 x10~ 5 


3.123x10~ 4 


COv(D A{l BB,D AA Bb) 


3.278 x10~ 4 ±7.803x10~ 5 


3.217x10~ 4 


CO-v(D Aa BB,D AaB b) 


-3.310x10~ 4 ±9.101 x10~ 5 


-3.366 x10~ 4 


cov(D AAB b,D AaB b) 


-4.654 x 1 0~ 4 ±5.826 x 1 0~ 5 


-5.191 x10~ 4 



*Parameter settings are the immigration rate of pollen m P = 0.08 and seeds m 5 = 0.04, the effective population size = 100, the selection coefficients in the gametophyte 
and sporophyte stages s aO = 5 o p = s a /2 = 0.04, and the degree of dominance =0.5. The genotypic frequencies in the continent and the initial island populations are 
0.1 225 for AABB, 0.1 05 for AABb, 0.0225 for AAbb, 0.1 05 for AaBB, 0.245 for AB/ab, 0.045 for ab/aB, 0.1 05 for Aabb, 0.0225 for aaBB, 0.1 05 for aaBb, and 0.1 225 for aabb. 
The steady-state simulation results (mean ± S d ) are obtained from 10000 independent runs. 
doi:1 0.1 371 /joumal.pone.0080538.t003 



selection coefficient (Figure 5c). The covariances between the 
zygotic LDs for the genotypes with a common genotype at the 
selective locus, i.e. E{cov{D AABB , Awsi)) and E(codD AaBB , D AaBb )), 
are less sensitive to selection than the covariances between the 
zygotic LDs for the genotypes with a common genotype at the 
neutral locus, i.e. EicoviD^s, D AaBB )) and E^oviD^i,, D AaBb )). 

The above results indicate that the gametic LD can have a 
similarly changing pattern to some zygotic LDs with the selection 
pressure. This provides the genetic basis of using zygotic LDs to 
describe genetic hitchhiking effects at the diploid level. Further- 
more, the covariances between gametic and zygotic LDs or 
between distinct zygotic LDs are informative to indicate genetic 
hitchhiking effects. 

Neutral Process. How zygotic LDs evolve in a purely neutral 
process forms another important issue to study the pattern of 
genotypic diversity along chromosomes since most molecular 
population evolution is governed by the neutral process. This also 
provides theoretical perception of using zygotic LD to reveal the 
structure of genomic diversity. Suppose that both loci are subject 
to the balance among the effects of genetic drift, recombination, 
and immigration. All items with selection coefficients are removed 
in the formulae in Appendix S2, but all the formulae in Appendix 
S3 remain unaltered. To assess the steady-state zygotic LDs and 
other covariances (Appendices S2 and S3), I need to calculate the 
following fifty-four expectations: E(p A ~'p B ~')(i, j— 0, 1, 2, 3,4; 
except i=/ = 4), E{jfc l £p 'D AB ){i,j= 0,1,2,3), Etff'fc 1 ' ^{i, 
7 = 0,1,2), E(p A - i p B - J D 3 AB )( hJ = Q,l), and E{D AB ). These expec- 
tations can be numerically calculated with Matkematka tool in 
different equations or different groups of equations. 

Letting g=p' A (i= 1, 2, 3, 4) in Eq. (3), I can obtain 



■ 4NmQ A 

4Nm + 1 — 1 



(11) 



The zth moment of allelic frequency is the same as that derived 
under a neutral process for individual loci since LD does not affect 
allelic frequency distribution. Let .F = (l +4A r m) _1 (an inbreeding 
coefficient in the island population) [57], which is equal to 
(E(p 2 A )-(E(p A )) 2 )/Q A (l-Q A ), 
differentiation coefficient F st in 
plants [58]. Eq. (11) 
, (i-W-Q A )F 



analogous to the population 
the classical island model for 
can be rewritten as 

E(p' A l ), the same as Wright's 



E(P' A )= ■ l+{i _ 2 )F 

expression except for plant species here ([11], p.450). Eq. (11) 
represents the steady-state moments of allelic frequency under the 
balance of migration-genetic drift, different from Robertson's 
results in a progressive inbreeding process [59]. Replacement of 
subscripts A with B in Eq. (1 1) yields E(p B ) (i = 1, 2, 3, 4). 

Substituting g=p A p B and g = D AB separately in Eq. (3) to yield 
two equations, I can obtain 



E(D AB ) = 



4Nm(\-r)D AB 
4Nr + (\-r)(4Nm+\y 



E(p a Pb) = QaQb + 



(\-r) 2 D AB 
4Nr + (\-r)(4Nm+\y 



(12) 



(13) 
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-a a 

m — I 



0.08 
0.06 - 
0.04 
0.02 
0 

-0.02 
-0.04 



0.01 0.02 0.03 

s aP ° rS aO° rS c/ 2 

DAB 9 DAABB 

DAaBB & DAABb 

- DAaBb 



0.04 




-0.0002 H 1 1 , 

0 0.01 0.02 0.03 0.04 

s aP ° rs aO ° rs J 2 

— • — cov(DAB, DAABB) — * cob(DAB, DAaBB) 

A cov(DAB,DAABb) 1 cov(DAB,DAaBb) 



0.0004 



0.0002 
0 

-0.0002 - 
-0.0004 - 
-0.0006 



0 



0.01 



0.02 



0.03 



0.04 



cov( DA ABB, DAaBB) 
cov(DA ABB, DAaBb) 
cov(DAaBB,DAaBb) 



s ao ors/2 

— 9 cov(DAABB.DAABb) 

1 cov(DAaBB, DAABb) 

cov(DAABb, DAaBb) 



Figure 5. Genetic hitchhiking effects on the steady-state 
gametic and zygotic LDs, and other types of covariances. 

Gametic and zygotic LDs (a); covariances between gametic and zygotic 
LDs (b); and covariances between distinct zygotic LDs (c). Results are 
obtained from the analytical model in the section of Analytical Theory. 
Parameter settings are the immigration rate of pollen m P = 0.08 and 
seeds m s = 0.04, the effective population size = 100, the selection 
coefficients in the gametophyte and sporophyte stages s a0 = s aP = s a l 
2, and the degree of dominance = 0.5. The genotypic frequencies in the 
continent and initial island populations are 0.1225 for AABB, 0.105 for 
AABb, 0.0225 for AAbb, 0.105 for AaBB, 0.245 for AB/ab, 0.045 for ab/aB, 
0.105 for Aabb, 0.0225 for aaBB, 0.105 for aaBb, and 0.1225 for aabb. 
doi:1 0.1 371 /journal.pone.0080538.g005 



Eq. (12) indicates that the expectation of gametic LD is equal to 
zero in the absence of LD in migrants, such as in a completely 
isolated population. Eq. (13) indicates that the expectation of joint 
allele frequencies at two loci is related to the gametic LD in 
migrants (Dab # 0) although expectations of individual allele 
frequencies are independent from each other (Eq. (1 1)). 

Expectations of the remaining forty-four functions can be 
calculated in the following steps. Substitution of g in Eq. (3) by 
functions Pa 'Pb ''^' j = 0, 1 ; except i=j=l), p x A l p X B 'D AB 

7 = 0,1; except i=j—l), and D AB can yield seven equations 
that can be used to numerically calculate their expectations: 
E(PA-'p 2 B~ J )(hJ = 0A; except i=j= 1), E( 1 ,\-p l B J D AB )( l ,j = Q,\- 
except i—j = 1), and E(D^ B ). Substitution of g in Eq. (3) by 
functions jfofc \i = 1 ,2), 'p\(i =1,2), faff 'D AB (i =1,2), 
p A ~'p B D AB (i=l,2), p A D 2 AB , and PbD 2 ab , yields ten equations to 
calculate their expectations. Substitution of g in Eq. (3) with two 
functions p\pb an d p A D AB yields two equations to calculate 
E(p\ps) and E(p A D AB ). Substitution of g in Eq. (3) with two 
functions p A p\, and p\D A B yields two equations to calculate 
E(PaP b ) and E(p B D AB ). Substitution of g in Eq. (3) with functions 
PaPb, PaPb, PapIDab, p\pbD ab , P A D 2 AB , and p 2 B D 2 AB , yields six 
equations to calculate their expectations: E(p 2 A p B ), and 
E(p 2 B D 2 AB ). Finally, substitution of g in Eq. (3) with the remaining 
seventeen functions p A p B , and D AB , yields seventeen equations 
to calculate their expectations: E(p A p B ), and E(D AB ). 

The above order of g substitutions with different functions is 
sequentially arranged since calculations of the expected functions 
in the later equations need the expectations of the functions 
derived from the former equations. All these calculations can be 
done using Mathematica equation solution tool. Expectations of 
high-order LD functions, such as E(D 2 AaBh ), can also be calculated 
with additional equations by setting g in Eq. (3) with different 
functions. These are not explored further. 

Once the expectations of the above fifty-four functions are 
available, the expectations of lower or the same order functions 
can be indirectly calculated. For instance, I can obtain 
E(p AB ) = E(p A p B ) + (l-r)E(D AB ) and E(p* AB p* ab ) = E(p A p a p B p h ) 
+ (l-r){E(p A p B D AB ) + E(j> a p h D AB )) + (l-r) 2 E(D AB ). The ex- 
pectations of all possible zygotic LDs, the variances of zygotic LDs, 
and the covariances among different LDs at the steady state can be 
calculated according to Appendices S2 and S3. 

For instance, the expectation of steady-state zygotic LDs for the 
genotype with double heterozygotes (E(D AaB b) from 
E(M(S DAaBh )) = 0), its variance (E{V(D AaBb )) = E(V(8 DAM ))), 
and its covariance with gametic LD (E(cov(D AB ,D AuB i,)) 
= E(cov(S Dab ,3 Dam ))) are given by 

E(D AaBh ) = 2(\ -r)E((l -2p A )(\ -2p B )D AB ) + 4(1 -r) 2 E(D 2 AB ) 
+ 2(m s + rap) (4E(p A p a p B p b ) - E(p* AB p* ah ) - E(p* Ab p* aB )) 

+^QABE(p: h )+Q: h E(p AB )+Q Ah E( P : B )+Q: B E(p Ah ) 

- 2{Q A E{p a p B p b ) + Q a E(p A p B Pb) 
+ Q,BE(p A p a pb) + QbE{p A p a p B )) 
+ m s (Q Aa Bb - 2Q Aa E(p B pb ) -2Q B bE(p A p a )) 



(14) 
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Table 4. Comparison between the simulation results and analytical model predictions under a neutral process* 





MC simulations 


Analytical model 


PA 


0.5079±0.0898 


0.5 


Pb 


0.5102±0.0913 


0.5 


Dab 


0.0974±0.0399 


0.0615 


Daabb 


0.0623 ±0.0363 


0.0360 


DaiiBB 


-0.0251 ±0.0312 


-0.0100 


D AaABb 


-0.0256±0.0312 


-0.0100 


DauBH 


0.0462 ±0.0383 


0.0200 


cov(D A b,Daabb) 


1.671 x10~ 4 ±3.842x10~ 5 


1.689x10~ 4 


co\(D A b,D Aii bb) 


-0.909 x10~ 4 ±4.71 x10~ 5 


-0.634 x10~ 4 


COv(D AB ,D A ABb) 


-0.916 x10~ 4 ±4.67x10~ 5 


-0.634x10~ 4 


CO\(D A B,D A uBb) 


1.774x10~ 4 ±5.876x10~ 5 


1.269x10~ 4 


CO\(D A ABB,D A aBB) 


-3.164x10~ 4 ±8.517x10~ 5 


-2.921 x10~ 4 


cov(D A abb,£>aabi,) 


-3.167x10~ 4 ±8.527x10~ 5 


-2.921 x10~ 4 


CO\(D A ABB,D A aBb) 


2.318x10~ 4 ±8.527x10~ 5 


2.193x10~ 4 


COv(D A{l BB,D A ABb) 


2.575x10~ 4 ±8.78x10~ 5 


2.267 x10~ 4 


CO-v(D Aa BB,D Aa Bb) 


-2.977 x10~ 4 ± 1.002 x10~ 4 


-3.237 x10~ 4 


mv(D A ABb,D A aBb) 


-2.949 x 1 0~ 4 ±0.990 x 1 0~ 4 


-3.237 x10~ 4 



*Parameter settings are the immigration rate of pollen m P =0.08 and seeds m 5 = 0.04, the effective population size = 100, and the recombination rate =0.05. The 
genotypic frequencies in the continent and the initial island populations are 0.1225 for AABB, 0.105 for A ABb, 0.0225 for AAbb, 0.105 for AaBB, 0.245 for4fi/af>, 0.045 for 
ab/aB, 0.105 for Aabb, 0.0225 for aaBB, 0.105 for aaBb, and 0.1225 for aabb. The steady-state simulation results (mean ± S^) are obtained from 10000 independent runs. 
doi:1 0.1 371 /journal.pone.0080538.t004 



E( V(D AaBh )) = — (4E(p A p a p B p h ( 1 - 2p A p a )( 1 - 2p B p h )) 

+ 2(1 -r)E((l -2 PA ){\ -4p APa )(l -2p B )(l -Ap B p b )D AB ) ( 15 ) 

- 16(1 -rfE{(\ -2p A )(\ -2p B )D AB )-\6(\-rfE{D\ B )) 

E(Cov(D AB ,D AaBh )) = ~ (E(p A p a p B p h (\ - 2p A )(\ - 2p B )) 

+ (1 -r)E({\ -2p A f(\ -2p A ) 2 D AB ) (16) 
+ 4(1 -r)E(p A p a p B p h D AB ) 

- ( 1 - r) 2 E({\ - 2p A )( 1 - 2p B )D 2 AB ) -4(1 - rfE{D AB )) 

Simulations confirm that the above analytical model performs 
well. The gametic and zygotic LDs and other covariances between 
two neutral loci can quickly reach steady-state distributions, 
reflecting the equilibrium among the effects of migration, 
recombination, and genetic drift. All analytical results are 
distributed within the range of one standard deviation of the 
simulation results (Table 4). Simulations also confirm that the 
expectations of Z>41S4 an d D AaBB and their covariances with 
gametic or other zygotic LDs are the same because both loci are 
neutral. This symmetry may help to reduce the number of 
expectations of distinct functions in theory. 

Figure 6 shows that different patterns exist for the expectations 
of gametic and zygotic LDs and other covariances with the 
recombination rate. E(D AB ) decreases faster than the absolute 
expectations of zygotic LDs with the recombination rate in 
addition to their inequality in magnitude (Figure 6a). Figure 6b 
shows that the absolute expectations of the covarainces between 



gametic and zygotic LDs gradually decrease with the recombina- 
tion rate for the genotypes with heterozygotes at one locus or two 
loci. EicoriDjus, A4ABb)) slightly decreases with the recombination 
rate, but does not approach zero in the presence of immigration 
that maintains gametic and zygotic LDs. Figure 6c shows that the 
covariances between different zygotic LDs are generally not as 
sensitive as some covariances between gametic and zygotic LDs to 
the change of linkage distance. E(cov(D AABB , E> AaBb )) and E(cov 
{D AaBB , AlAti*)) slightly decrease with the recombination rate, while 
the covariances between the zygotic LDs of the genotypes with one 
common genotype at a locus slightly increase with the recombi- 
nation rate. 

The above results indicate that a neutral process can generate a 
similar pattern between zygotic and gametic LDs along chromo- 
somes, with strong LDs within short distances and weak LDs 
within long distances. The covariances between gametic and 
zygotic LDs or between distinct zygotic LDs are relatively 
insensitive to the linkage distance. 

Discussion 

In this study, I have developed the evolutionary theory of 
zygotic LDs in a local plant population, complementing the 
previous theories that mainly focus on the statistical issues [5] , [6] , 
[9]. The theory shows that evolutionary forces can generate 
different patterns among gametic and zygotic LDs, the covariances 
between gametic and zygotic LDs, and the covariances between 
different zygotic LDs. Zygotic LDs can be greater or smaller than, 
or comparable to gametic LD, depending on the major ecological 
and evolutionary processes involved in a local population. The 
covariances between gametic and zygotic LDs are more sensitive 
to the effects of mating system, linkage distance, and genetic drift, 
than to the effects of seed and pollen flow and selection. The 
covariances between different zygotic LDs are relatively robust to 
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Figure 7. A comparison of transient gametic and zygotic LDs in 
a finite isolated population versus in an infinite population. 

Zygotic LDs in the finite population are calculated by synthesizing the 
theories of Robertson [59] and Ohta and Kimura [12]: N= 10, f = 10, and 
N=10, f = 20 in (a). Gametic LD in the finite population is calculated 
from Hill and Robertson [8]: W=10, f=10, and W=10, f = 20 in (b). 
Gametic and zygotic LDs in the infinite population are calculated from 
Weir and Cockerham [4]: a = 0.05, f=10; a = 0.05, f = 20; a = 0.95, f=10; 
and a = 0.95, f = 20 in (a) and (b). The initial settings for the finite 
population are W=10 and the frequency of double heterozygotes 
(coupling) =1 (gametic LD = 0.25 and the allelic frequency at each 
diallelic locus = 0.5). The initial setting for the infinite population is the 
frequency of double heterozygotes (coupling) =1 (gametic LD = 0.25 
and the allelic frequency at each diallelic locus = 0.5). 
doi:10.1371/journal.pone.0080538.g007 



Figure 6. Effects of the linkage distance on the steady-state 
gametic and zygotic LDs and other types of covariances in a 
neutral process. Gametic and zygotic LDs (a); covariances between 
gametic and zygotic LDs (b); and covariances between distinct zygotic 
LDs (c). Results are obtained from the analytical model in the section of 
Analytical Theory. Parameter settings are the immigration rate of pollen 
m P = 0.08 and seeds m s = 0.04, the effective population size=100, the 
selection coefficients in the gametophyte and sporophyte stages 
s oo = s ap = s a = 0, and the degree of dominance =0.0. The genotypic 
frequencies in the continent and initial island populations are 0.1 225 for 
AABB, 0.105 for AABb, 0.0225 for AAbb, 0.105 for AaBB, 0.245 for AB/ab, 
0.045 for ab/aB, 0.105 for Aabb, 0.0225 for aaBB, 0.105 for aaBb, and 
0.1225 for aabb. 

doi:1 0.1 371 /journal.pone.0080538.g006 



the change in gene flow, selection, and genetic distance, but 
sensitive to the genetic drift effects and mating system. Consistent 
patterns exist for the covariances between zygotic LDs for the 
genotypes with a common genotype at one locus, or for the 
genotypes without any common genotype at each locus. These 
similarities and differences suggest the potential utility of zygotic 
LDs in revealing the ecological and evolutionary processes 
underlying the pattern of population genomic diversity at the 
diploid level. 

It is important to understand that in a pure drift process, LD is 
transient in a completely isolated population of random mating. 
Expectations of both gametic and zygotic LDs are zero although 
the expectations of their squared values are nonzero [8], [44], 
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[60]. Zygotic LDs are smaller than gametic LDs in magnitude, but 
decay more slowly than gametic LDs with time [12], [46]. This is 
because the genotypic association may primarily arise from the 
effects other than the linkage distance although the linkage 
distance can affect the frequencies of genotypic recombinants. 
Weir and Cockerham[4] and Cockerham and Weir [3] have 
decomposed gametic and zygotic LDs in terms of descent 
measures for an infinite population with a mixed mating system. 
Both gametic and zygotic LDs decay faster in an infinite than in a 
finite population within a short linkage distance when the genetic 
drift effects are in the same order as the selfing rate (al/2N e ; e.g., 
Figure 7). A predominandy selfing population reduces the rates of 
decay of both gametic and zygotic LDs. When additional driving 
forces are involved, the above "null" expectation and the rates of 
decay in gametic and zygotic LDs could be changed [1],[2]. 

Note that the theory only addresses the constant immigration of 
seeds and pollen. In reality, a frequent situation is the stochastic 
migration of seeds and pollen due to the influences of biotic and 
abiotic factors [61], [62]. This occurs particularly when the source 
populations or the pollen and seed pools are unstable. Under this 
situation, the gametic and genotypic frequencies fluctuate in 
migrating seeds and pollen, and so do the gametic and zygotic LDs 
in migrants. Zygotic and gametic LDs can exhibit more 
fluctuations under the joint effects of genetic drift and stochastic 
migration. This can weaken the relationships between zygotic and 
gametic LDs. Nevertheless, the explored qualitative relationships 
between zygotic LDs and migration remain valid. How the 
stochastic migration of seeds and pollen affects the relative gametic 
and zygotic LDs remains unclear, and this forms a topic for further 
study. 

Also, note that a plant mating system in a natural population 
may exhibit a dynamic property [25]. Mating system can be 
naturally changed through different ways [63], such as the change 
of pollen pool and the shift from wind to animal pollinations [64] , 
[65]. Since zygotic LDs and other covariances are sensitive to the 
change of mating system (Figure 1), an unstable mating system 
enhances the fluctuation of these covariances. Nevertheless, the 
non-linear relationships between the selfing rate and zygotic LDs 
remain valid. 

Apart from the above assumptions, the theory suggests several 
useful implications [1], [2]. First, the newly explored relationships 
between gametic and zygotic LDs under different evolutionary 
forces, not the purely statistical relationships [3], [5], [6], suggest 
their different or similar practical applications. Different patterns 
between gametic and zygotic LDs indicate that they can be 
applied for different purposes. Zygotic LDs provide additional 
information for inferring population history. Previous studies 
emphasize the use of gametic LD for this purpose [61], [66]. The 
present theory shows that zygotic LDs exhibit more diverse 
patterns in response to different driving forces, which can reinforce 
our inference on the major ecological and evolutionary processes. 

The occurrence of a weak gametic LD combined with strong 
zygotic LDs suggests epistatic interactions at the diploid level (e.g., 
postzygotic isolation due to the Dobzhansky-Muller incompatibil- 
ity [27], [28], [47]). The occurrence of a strong gametic LD 
combined with weak zygotic LDs suggests the involvement of non 
epistatic processes at the diploid level, including migration, linear- 
additive selection, and genetic drift processes. When both strong 
gametic and zygotic LDs arise from the tight linkage, they can be 
applied for the same purposes. For example, in analyzing 
normalized gametic and zygotic LDs in a human population, 
the same SNP markers exist when both gametic and zygotic LDs 
are very strong (say, the squares of normalized gametic and zygotic 
LDs >0.9; [67]). The relatively stable patterns in zygotic and 



gametic LDs and in covariances between gametic and zygotic LDs 
across multiple populations suggest the impacts of seed and pollen 
flow or weak selection. Patterns from multiple samples of a given 
population or from multiple different natural populations can 
strengthen such inferences. 

Second, the theory provides a genetic basis of using zygotic LDs 
for QTL mapping that has been recendy addressed [23]. A similar 
pattern between zygotic and gametic LDs with the linkage distance 
implies the common utility for QTL mapping. Zygotic LD-based 
QTL mapping can be conducted in nonrandom mating popula- 
tions [23]. One caution is that spurious and unstable non-random 
associations can occur in natural populations under the influences 
of the driving forces other than the recombination process. This 
can influence the accuracy and precision of QTL mapping. QTL 
mapping based on the linkage maps from a single family, such as a 
half-sib family from a single tree or a full-sib family from a single 
cross, is not affected. However, the population-based linkage maps 
could be affected although this approach is commonly suggested to 
search for LDs within a short linkage distance at a finer scale [68], 
[69]. Thus, the patterns of zygotic LDs can be used to 
preliminarily screen markers for QTL mapping through a high 
criterion [46], or to effectively remove spurious LDs through a 
deliberate experiment [70] . This may improve QTL mapping with 
the population-based linkage maps. 

Third, the theory aids in predicting the effects of seed and pollen 
flow on zygotic LDs in a local population. Previous studies use 
gametic LD to estimate gene flow in a specific case, such as in 
hybrid zones [71]. The present theory shows that gametic LD is 
more sensitive than zygotic LDs to either seed flow or pollen flow. 
Seed flow has greater effects than pollen flow on gametic LD. In 
natural populations of flowering plants, pollen flow is often more 
extensive than seed flow among mature populations, especially for 
the predominandy outcrossing species [72]. The cumulative effects 
on gametic LD from pollen flow could be substantial. The robust 
pattern of zygotic LDs to the impacts of seed or pollen flow enables 
their utility for inferring if gametic LD is generated by the forces 
other than migration. One extreme case is the admixture of two or 
more plant or animal populations, such as cross breeding, which 
results in the same consequence as that produced by a large 
proportion of immigrating seeds. This produces extensive gametic 
LDs rather than zygotic LDs [61], [66]. Only those tighdy linked 
loci can maintain strong zygotic and gametic LDs [46] . Thus, the 
multilocus patterns of joint gametic and zygotic LDs can be used 
to judge if immigration is an important process to shape gametic 
LDs in local populations. 

Fourth, the theory aids in assessing the selection mode (additive 
or epistatic) in the gametophyte and sporophyte stages in 
generating gametic and zygotic LDs. "Bulmer effects" mainly 
emphasize the impacts of selection on gametic LD [26], but 
gametic LD does not provide the information on the genotypic 
interaction at the diploid level. Extensive reports are recorded in 
the literature about the use of gametic LD for detecting selection 
signature along chromosomes [73]. So far, zygotic LDs have not 
been applied to detecting the genetic basis of adaption at the 
diploid level. In the linear additive-viability model, selection from 
the two stages is compounded. Gametic LD is greater than zygotic 
LDs in magnitude because selection affects gametic LD at each 
stage but affects zygotic LDs only in the sporophyte stage, similar 
to the effects of haploid pollen and diploid seed flow. However, in 
the presence of epistatic selection at the diploid level, some 
genotypes have zygotic LDs larger than gametic LD while other 
genotypes have zygotic LDs smaller than gametic LD. Such 
divergent patterns can aid in our inference on epistatic selection. 
One typical situation is a natural hybrid zone (a tension zone) [29] 
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where epistatic selection can cause zygotic LDs greater than 
gametic LD [41], which provides the information complementary 
to two non-allele interaction at the haploid level [27], [28], [30]. 
The joint patterns of gametic and zygotic LDs can be used to infer 
the selection mode (additive or epistatic) at the diploid level. 

In addition, the genotypic interaction on fitness may arise from 
the dominance by dominance effects for D AaBb , or the additive by 
dominance effects for T>AABb or D AaBB , or the additive by additive 
effects for D AABB at two loci. One further study is to assess the 
genetic mechanisms of these epistases in distinct zygotic LDs at the 
sporophyte stage. 

Finally, it is of interest to discuss the utility of the covariances 
between distinct zygotic LDs since few studies have examined such 
high-order LDs [9], [10]. The present theory suggests one robust 
property of these high-order LDs, i.e. the presence of a consistent 
pattern for the genotypes with one common genotype at one locus 
or for the genotypes without any common genotype at each locus. 
This property can be used to effectively determine the impacts 
from migration, recombination, and additive weak selection, and 
to assess the effects of effective population size and/ or a mating 
system. Given a stable effective population size and a stable mating 
system, a significant bias from the robust property implies epistatic 
selection (Table 2) or very diverse selection systems among 
genotypes. This requires further empirical verification with 
appropriate data collections. 
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